/*
 * 蒙特卡洛积分：求单位半圆面积
 */

#include<iostream>
using namespace std;
#include<cmath>

int main(){

	int i, n=1E8;
	double x, f, sum_f = 0, sum_f2 = 0, ave_f, ave_f2, I1, I2;

	for(int i=0;i<n;i++){
		x = (double)rand()/RAND_MAX*2 - 1;//取随机数 x 在 [-1,1] 均匀分布
		f = sqrt(1-x*x);// f*f + x*x =1
		sum_f += f;
		sum_f2 += f*f;
	}
	ave_f = sum_f / n;
	ave_f2 = sum_f2 / n;

	I1 = 2 * ( ave_f - 4 * sqrt( (ave_f2 - ave_f * ave_f)/n ) );//置信区间左边界
	I2 = 2 * ( ave_f + 4 * sqrt( (ave_f2 - ave_f * ave_f)/n ) );//置信区间右边界
	cout<<" If is in ["<<I1<<","<<I2<<"] by "<<99.993666<<"% chance."<<endl;
	// 4 sigma 置信度为 99.993666%
	cout<<" Exact value of If = "<<M_PI/2<<endl;

	return 0;
}
